#include <math.h>
#include <stdlib.h>
#include <time.h>

#include <chrono>
#include <iostream>

#include "Interpolation1D.h"
#include "Interpolation2D.h"
#include "Interpolation3D.h"
#ifdef _USEGSL
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#endif

typedef std::chrono::high_resolution_clock Clock;

#define RESET "\033[0m"
#define BLACK "\033[30m"  /* Black */
#define RED "\033[31m"    /* Red */
#define GREEN "\033[32m"  /* Green */
#define YELLOW "\033[33m" /* Yellow */
#define BLUE "\033[34m"   /* Blue */
#define PURPLE "\033[35m" /* Purple */
#define CYAN "\033[36m"   /* Cyan */
#define WHITE "\033[37m"  /* White */

#define NORMAL (std::cout << RESET)
#define INFO (std::cout << GREEN)
#define WARN (std::cout << YELLOW)
#define ERROR (std::cout << RED)
#define DEBUG (std::cout << CYAN)
#define END (std::endl)
#define REND "\033[0m" << std::endl

#ifdef _USEGSL
void test_gsl_1d(const gsl_interp_type *T, int N, int M) {
  srand(time(0));
  const int m = 7, n = 7;
  auto t1 = Clock::now();
  auto t2 = Clock::now();
  const double xa[] = {0.00117378, 0.00472397, 0.0171859, 0.0409838,
                       0.0780596,  0.130081, 0.198535};
  const double ya[] = {0.365616,  1.893040,  3.786090, 7.572170,
                       15.144300, 30.288700, 60.577400};

  double t_construct = 0.0, t_evaluate = 0.0;
  for (int i = 0; i < N; i++) {
    srand(time(0));
    t1 = Clock::now();
    gsl_spline *spline = gsl_spline_alloc(T, m);
    gsl_spline_init(spline, xa, ya, m);
    gsl_interp_accel *acc = gsl_interp_accel_alloc();
    t2 = Clock::now();
    t_construct +=
        std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();

    t1 = Clock::now();
    float x, y = 0.0;
    for (int j = 0; j < M; j++) {
      x = (((unsigned int) (rand() * y)) % 100 + 2) / 1000.0;
      y = gsl_spline_eval(spline, x, acc);
    }
    t2 = Clock::now();
    t_evaluate +=
        std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
    gsl_spline_free(spline);
    gsl_interp_accel_free(acc);
  }

  std::cout << "average construct time (for N = " << N << " times): " << GREEN
            << t_construct * 1.0 / N << RESET << " microseconds" << RESET
            << std::endl;
  std::cout << "average evaluate time (for N*M = " << N * M
            << " times): " << GREEN << t_evaluate * 1.0 / N / M << RESET
            << " microseconds" << RESET << std::endl;
}

void test_gsl_2d(const gsl_interp2d_type *T, int N, int M) {
  srand(time(0));
  const int m = 7, n = 7;
  auto t1 = Clock::now();
  auto t2 = Clock::now();
  const double xa[] = {0.00117378, 0.00472397, 0.0171859, 0.0409838,
                       0.0780596,  0.130081,   0.198535};
  const double ya[] = {0.365616,  1.893040,  3.786090, 7.572170,
                       15.144300, 30.288700, 60.577400};
  const double za[] = {0.0217822, 0.0253224, 0.0288237, 0.0346827, 0.0448323,
                       0.0636086, 0.100366,  0.0233179, 0.0268545, 0.0303556,
                       0.0362159, 0.0463659, 0.0651426, 0.101902,  0.0296429,
                       0.0331470, 0.0366371, 0.0425000, 0.0526603, 0.0714467,
                       0.108208,  0.0402311, 0.0440292, 0.0477457, 0.0538394,
                       0.0641187, 0.0829203, 0.119654,  0.0511250, 0.0554077,
                       0.0595859, 0.0662932, 0.0771901, 0.0963434, 0.133061,
                       0.0625876, 0.0673198, 0.0719785, 0.0794046, 0.0910973,
                       0.110757,  0.147656,  0.0748282, 0.0800098, 0.0851434,
                       0.0933663, 0.106111,  0.126669,  0.163872};

  double t_construct = 0.0, t_evaluate = 0.0;
  for (int i = 0; i < N; i++) {
    srand(time(0));
    t1 = Clock::now();
    gsl_spline2d *spline = gsl_spline2d_alloc(T, m, n);
    gsl_spline2d_init(spline, xa, ya, za, m, n);
    gsl_interp_accel *xacc = gsl_interp_accel_alloc();
    gsl_interp_accel *yacc = gsl_interp_accel_alloc();
    t2 = Clock::now();
    t_construct +=
        std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();

    t1 = Clock::now();
    float x, y, z;
    for (int j = 0; j < M; j++) {
      x = (((unsigned int) (rand() * z)) % 100 + 2) / 1000.0;
      y = (((unsigned int) (rand() * z)) % 100 + 1) / 10.0;
      z = gsl_spline2d_eval_extrap(spline, x, y, xacc, yacc);
    }
    t2 = Clock::now();
    t_evaluate +=
        std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
    gsl_spline2d_free(spline);
  }

  std::cout << "average construct time (for N = " << N << " times): " << GREEN
            << t_construct * 1.0 / N << RESET << " microseconds" << RESET
            << std::endl;
  std::cout << "average evaluate time (for N*M = " << N * M
            << " times): " << GREEN << t_evaluate * 1.0 / N / M << RESET
            << " microseconds" << RESET << std::endl;
}
#endif

void test1D(Interpolation1DBase<double, float> &p1, int N, int M) {
  srand(time(0));
  int m = 7, n = 7;
  auto t1 = Clock::now();
  auto t2 = Clock::now();
  float xa[] = {0.00117378, 0.00472397, 0.0171859, 0.0409838,
                0.0780596,  0.130081,   0.198535};
  float ya[] = {0.365616,  1.893040,  3.786090, 7.572170,
                15.144300, 30.288700, 60.577400};
  double t_construct = 0.0, t_evaluate = 0.0;
  for (int i = 0; i < N; i++) {
    srand(time(0));
    t1 = Clock::now();
    p1.construct(xa, ya);
    t2 = Clock::now();
    t_construct +=
        std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();

    t1 = Clock::now();
    float x, y;
    for (int j = 0; j < M; j++) {
      x = (((unsigned int) (rand() * y)) % 100 + 2) / 1000.0;
      y = p1.evaluate(x);
    }
    t2 = Clock::now();
    t_evaluate +=
        std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
  }

  std::cout << "average construct time (for N = " << N << " times): " << GREEN
            << t_construct * 1.0 / N << RESET << " microseconds" << RESET
            << std::endl;
  std::cout << "average evaluate time (for N*M = " << N * M
            << " times): " << GREEN << t_evaluate * 1.0 / N / M << RESET
            << " microseconds" << RESET << std::endl;
}


void test2D(Interpolation2DBase<double, float> &p1, int N, int M) {
  srand(time(0));
  int m = 7, n = 7;
  auto t1 = Clock::now();
  auto t2 = Clock::now();
  float xa[] = {0.00117378, 0.00472397, 0.0171859, 0.0409838,
                0.0780596,  0.130081,   0.198535};
  float ya[] = {0.365616,  1.893040,  3.786090, 7.572170,
                15.144300, 30.288700, 60.577400};
  float za[7][7] = {0.0217822, 0.0253224, 0.0288237, 0.0346827, 0.0448323,
                    0.0636086, 0.100366,  0.0233179, 0.0268545, 0.0303556,
                    0.0362159, 0.0463659, 0.0651426, 0.101902,  0.0296429,
                    0.0331470, 0.0366371, 0.0425000, 0.0526603, 0.0714467,
                    0.108208,  0.0402311, 0.0440292, 0.0477457, 0.0538394,
                    0.0641187, 0.0829203, 0.119654,  0.0511250, 0.0554077,
                    0.0595859, 0.0662932, 0.0771901, 0.0963434, 0.133061,
                    0.0625876, 0.0673198, 0.0719785, 0.0794046, 0.0910973,
                    0.110757,  0.147656,  0.0748282, 0.0800098, 0.0851434,
                    0.0933663, 0.106111,  0.126669,  0.163872};
  float *zp[7];
  for (int i = 0; i < 7; i++) zp[i] = &za[i][0];

  double t_construct = 0.0, t_evaluate = 0.0;
  for (int i = 0; i < N; i++) {
    srand(time(0));
    t1 = Clock::now();
    p1.construct(xa, ya, zp);
    t2 = Clock::now();
    t_construct +=
        std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();

    t1 = Clock::now();
    float x, y, z = 0.0;
    for (int j = 0; j < M; j++) {
      x = (((unsigned int) (rand() * z)) % 100 + 2) / 1000.0;
      y = (((unsigned int) (rand() * z)) % 100 + 1) / 10.0;
      z = p1.evaluate(x, y);
    }
    t2 = Clock::now();
    t_evaluate +=
        std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
  }

  std::cout << "average construct time (for N = " << N << " times): " << GREEN
            << t_construct * 1.0 / N << RESET << " microseconds" << RESET
            << std::endl;
  std::cout << "average evaluate time (for N*M = " << N * M
            << " times): " << GREEN << t_evaluate * 1.0 / N / M << RESET
            << " microseconds" << RESET << std::endl;
}

void test3D(Interpolation3DBase<double, float> &p1, int N, int M) {
  srand(time(0));
  int m = 7, n = 7, l = 7;
  auto t1 = Clock::now();
  auto t2 = Clock::now();
  float xa[] = {0.00117378, 0.00472397, 0.0171859, 0.0409838,
                0.0780596,  0.130081, 0.198535};
  float ya[] = {0.065616,  1.893040,  3.786090, 7.572170,
                15.144300, 30.288700, 60.577400};
  float za[] = {0.003067, 0.006772, 0.050856, 0.165419, 0.269267, 0.431288, 1.07437};
  float ***va = new float **[m];
  for (int i = 0; i < m; i++)
  {
    va[i] = new float*[n];
    for (int j = 0; j < n; j++)
    {
      va[i][j] = new float[l];
      for (int k = 0; k < l; k++)
        va[i][j][k] = ((int) (rand())) % 100 / 10.0;
    }
  }

  double t_construct = 0.0, t_evaluate = 0.0;
  for (int i = 0; i < N; i++) {
    srand(time(0));
    t1 = Clock::now();
    p1.construct(xa, ya, za, va);
    t2 = Clock::now();
    t_construct +=
        std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();

    t1 = Clock::now();
    float x, y, z, v = 0.0;
    for (int j = 0; j < M; j++) {
      x = (((unsigned int) (rand() * v)) % 100 + 2) / 1000.0;
      y = (((unsigned int) (rand() * v)) % 100 + 1) / 10.0;
      z = (((unsigned int) (rand() * v)) % 100 + 1) / 100.0;
      v = p1.evaluate(x, y, z);
    }
    t2 = Clock::now();
    t_evaluate +=
        std::chrono::duration_cast<std::chrono::microseconds>(t2 - t1).count();
  }

  std::cout << "average construct time (for N = " << N << " times): " << GREEN
            << t_construct * 1.0 / N << RESET << " microseconds" << RESET
            << std::endl;
  std::cout << "average evaluate time (for N*M = " << N * M
            << " times): " << GREEN << t_evaluate * 1.0 / N / M << RESET
            << " microseconds" << RESET << std::endl;
}

int main(int argc, const char *argv[]) {
  int m = 7, n = 7, l = 7, N = 1000, M = 1000;
  std::cout << "-----------Time Test-----------" << std::endl;
  if (argc != 3) {
    std::cerr << "No N and M specified. Default N = 1000, M = 1000."
              << std::endl;
  } else {
    N = atoi(argv[1]);
	  M = atoi(argv[2]);
  }

  std::cout << "Constructing the Interpolation Class " << CYAN
            << "for N = " << N << " times. " << RESET << std::endl;
  std::cout << "For each construct, evalute the value " << CYAN
            << "for M = " << M << " times. " << RESET << std::endl;

  std::cout << CYAN << "1D Linear Interpolation" << RESET << std::endl;
  if (N * M != 0) {
    LinInt1D<double, float> l1 = LinInt1D<double, float>(m);
    test1D(l1, N, M);
  }
#ifdef _USEGSL
  const gsl_interp_type *T1 = gsl_interp_linear;
  std::cout << PURPLE << "GSL 1D linear Interpolation" << RESET << std::endl;
  test_gsl_1d(T1, N, M);
#endif

  std::cout << CYAN << "2D Linear Interpolation" << RESET << std::endl;
  if (N * M != 0) {
    LinInt2D<double, float> l2 = LinInt2D<double, float>(m, n);
    test2D(l2, N, M);
  }
#ifdef _USEGSL
  const gsl_interp2d_type *T2d = gsl_interp2d_bilinear;
  //const gsl_interp2d_type *T2 = gsl_interp2d_bicubic;
  std::cout << PURPLE << "GSL 2D linear Interpolation" << RESET << std::endl;
  test_gsl_2d(T2d, N, M);
#endif

  std::cout << CYAN << "3D Linear Interpolation" << RESET << std::endl;
  if (N * M != 0) {
    LinInt3D<double, float> l3 = LinInt3D<double, float>(m, n, l);
    test3D(l3, N, M);
  }

  std::cout << CYAN << "1D Quadratic Interpolation" << RESET << std::endl;
  // std::cin >> N >> M;
  if (N * M != 0) {
    QuaInt1D<double, float> p1 = QuaInt1D<double, float>(m);
    test1D(p1, N, M);
  }

  std::cout << CYAN << "2D Quadratic Interpolation" << RESET << std::endl;
  // std::cin >> N >> M;
  if (N * M != 0) {
    QuaInt2D<double, float> p2 = QuaInt2D<double, float>(m, n);
    test2D(p2, N, M);
  }

  std::cout << CYAN << "3D Quadratic Interpolation" << RESET << std::endl;
  // std::cin >> N >> M;
  if (N * M != 0) {
    QuaInt3D<double, float> p3 = QuaInt3D<double, float>(m, n, l);
    test3D(p3, N, M);
  }

  std::cout << CYAN << "1D Spline Interpolation" << RESET << std::endl;
  // std::cin >> N >> M;
  if (N * M != 0) {
    SplInt1D<double, float> s1 = SplInt1D<double, float>(m);
    test1D(s1, N, M);
  }
#ifdef _USEGSL
  const gsl_interp_type *T2 = gsl_interp_cspline;
  std::cout << PURPLE << "GSL 1D spline Interpolation" << RESET << std::endl;
  test_gsl_1d(T2, N, M);
#endif

  std::cout << CYAN << "2D Spline Interpolation" << RESET << std::endl;
  // std::cin >> N >> M;
  if (N * M != 0) {
    SplInt2D<double, float> s2 = SplInt2D<double, float>(m, n);
    test2D(s2, N, M);
  }

  std::cout << CYAN << "3D Spline Interpolation" << RESET << std::endl;
  // std::cin >> N >> M;
  if (N * M != 0) {
    SplInt3D<double, float> s3 = SplInt3D<double, float>(m, n, l);
    test3D(s3, N, M);
  }
}
